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1 Abstract: Nonlinear elliptic system for generating adaptive quadrilateral meshes in curved domains 

■ is presented. Presented technique has been implemented in the C ++ language. The included software 

package can write the converged meshes in the GMV and Matlab formats. Since, grid adaptation is 
required for numerically capturing important characteristics of a process such as boundary layers. So, 
the presented technique and the software package can be a useful tool. 
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g 1 Introduction 

Quadrilateral grids are extensively used for numerical simulation. Accuracy of a simulation is 
j> \ strongly depend on the grid quality. Here, quality means orthogonality at the boundaries and 

quasi-orthogonality within the critical regions, smoothness, bounded aspect ratios and solution 
adaptive behaviour. Grid adaptation is used for increasing the efficiency of numerical schemes 
by focusing the computational effort where it is needed. In this article, we review the elliptic 
grid generation system for generating adaptive quadrilateral meshes. The presented scheme has 
been implemented in the C ++ language. 

For meshing a domain into non- simplex elements (quadrilaterals in 2D and hexahedrals 
in 3D), we seek a mapping from a reference square or cube to the physical domain. This 
mapping can be algebraic in nature such as Transfinite Interpolation or it can be expressed 
by a system of nonlinear partial differential equations JH; 0, and references therein] such as 
elliptic system. We are looking for a vector mapping, = (^,y) ? , from a unit square in 

the reference space (k = [0, 1] x [0, 1]) to a physical space (k); i.e. k 1 — ► k (see Figured])- 
Mapping gives the position of a point in the physical space corresponding to a point in the 
computational or reference space. Let the physical space be given by the x and y coordinates 



and the computational space be given by the t, and r\ coordinates (<§ G [0, 1] and r\ E [0, 1]). We 
are using the following elliptic system for defining the mapping j^t = (x, y)' 



d 2 y d 2 y d 2 y 

822 W ~ 2gil W^i +8n d¥ +Pn+Qyn = °- 

Here, the terms P and Q are used for grid adaptation and are given as 



(1) 

(2) 



P = g22 Pn -2g 12 Pl 2 + gnPL (3) 
Q = S22 Pn -2g n P 2 2 + g n P 2 2 . (4) 

Equations [T]|21 are non-linear and are coupled through the metric coefficients gij (coefficients of 
the metric tensor). Metric coefficients are given as 

gn=*!+y|, g22=x 2 +y 2 and gn=x^x n +y^y T} . (5) 

For generating grids in the physical space, the elliptic system HE is solved for the coordinates 
(x,y) on a unit square in the computational space by the method of Finite Differences. Boundary 
of the physical domain is specified as the Dirichlet boundary condition on the unit square in the 
computational space. In the Figure [H gi (= r«) and g2 (= r^) are the covariant base vectors 
at the point (xi,yj). Figure El shows a finite difference stencil around the point (£,-, 7]j) in the 
computational space. A finite difference approximation of and x-q at the point (i, j) (see 
Figure |2j) is 

[x(i+l,j)-x(i-l,j)] [x(i,j+l)-x(i,j-l)] 

x t — r - TT an d Xn = — . 

Similarly, x^ and y^ can be defined. Here, we are assuming that the grid in the computational 
space is uniform. However, grid in the physical space can be compressed or stretched. Terms 
Pfj (i =1,2 and j = 1,2 and k = 1,2 and P\ 2 = P 2 i) in the equations O-© are determined 
through another mapping &\. The mapping &\ is shown in the Figure This mapping maps 
a unit square in the computational space to a unit square in the parameter space. For defining 
the mapping &\ : k — > k\ , the boundary and internal grid points of the parameter space are 
mapped to the reference space. The Jacobian matrix T of the mapping J^i and the vectors Pn, 
P12 and P22 are given as follows 

The terms Pjj (i, j = 1,2) are the first component of the vector P ;; - and the terms Pfi are the 
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Figure 1 : Mapping ^ k from a reference unit square (k) on the left to a physical domain (k). 
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Figure 2: Finite difference stencil in the ^-J] computational space. 



second component of the vector P z -j. It should be noted that the vectors Pn, P12 and P22 can 
be computed a priori for clustering the grid points in the physical space. A second order finite 
difference approximation of different operators required for computing the vectors Pn, P22, P12 
and the Jacobian T are given in the Table [TJ We are using the stencil shown in the Figure |2j 



2 C ++ Implementation 

We have implemented the presented technique in the C ++ language for generating adaptive 
grids. The package can write meshes in the Matlab and GMV 11311 formats. It consists of one 
Domain class (see the subsections l2.2l and l2.3l) . Domain class is used for expressing unit square 
in the computational space (see line no. 045 in the subsection l2.ll) . unit square in the parameter 
space (see line no. 022 in the subsection l2.ll ) and the physical domain (see line no. 038 in the 
subsection l2.ll ). The physical domain is defined in the file functions.h (see the sub section 12 .61) . 
For clustering grids in the parameter space different functions are defined in the domain class 
(see the line numbers 026, 027, 029, 030, 032, 033 in the subsection l2~2b. 

The coupled elliptic system are linearised by the method of Finite Difference and the re- 
sulting system is solved by the SOR relaxation (see the subsection 12.5b . The SOR algorithm 
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Figure 3: Mapping &\ from a unit square (k) in the reference space to a unit square in the 
parameter space (k\). 



Table 1: Finite difference approximation of continuous operators. 

s(i + lj)-s(i-l, j) s(i + 1J) - 2s(i, j) + s(i - 1 , j) 
S ^ = 2A| ' *« = 

_ + i)-t(i, J-i) . _ j + 1) - 2f(i, y) y- 1) 

2A^ ' *™- Arp 

_ 5(/,j + i)-25(/,j)+5(/,j-i) _ f (? + 1 , y) - 2 f j) + f - 1 , y) 
5,7,7 " a^ ' X? 

_ s(mj + l)+s(z-lj-l)-s(/-l,j + l)-s(/+l,,/-l) 
S§ " 4A<^At] 
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consists of three loops, while outer loop (see line number 043 in the subsection 12.51) and two 
inner for loops (see line numbers 046 and 047 in the subsection l2.5l) . Each iteration of an inner 
loop provides a new mesh by the SOR relaxation. The outer loop is controlled by the maximum 
number of SOR iterations (see line number 028 in the sub section 12 .5 1) and a given tolerance (see 
the line number 027 in the subsection !2.5l) . 

The overall algorithm proceeds as follows. Generate grids in the computational and pa- 
rameter spaces. Compute the matrix T and vectors P;,- for defining the mapping &\ (from 
computational space to parameter space). An initial grid, y q \&, in the physical region is gener- 
ated (say by Transfinite Interpolation). This information is then passed to the SOR solver (see 
the line number 043 in the subsection l2.ll) . 



2.1 main.cpp 

ooi //+++++++++++++++++ 

02 # in elude <iost ream> 

003 # include <iomanip> 

004 #include <vector> 

005 #include <iterator> 
6 # include <f stream> 

007 #include <sstream> 

008 #include <map> 

9 # include "domain . h" 

010 // f include "write_matlab . h " 

011 # include "matrix . h" 

012 #include " sor_solver . epp" 

013 //+++++++++++++++++ 



014 


int 


main ( ) { 




015 




bool grid_dist = true ; 




016 




bool run_ellip = true; 




on 




unsigned xdim, ydim ; 




018 




xdim = 31, ydim = 31; 




019 




double del_xi = 1 . /double (xdim-1 . ) ; 




020 




double del_et a = 1.0/ double (ydim-1 . ) ; 




021 




//Parameter Space ref (xdim, ydim) ; 




022 




Domain parm {xdim, ydim) ; 




023 




/ /Meshing the Parameter 




024 




parm. Grid_Gen ( ) ; 




025 




//Clustering the Mesh 




026 




//Example 1 




027 




//parm.Cluster_X_Near (0.5) ; 




028 




//parm.Cluster_Y_Near (0.5) ; 




029 




//Example 2 




030 




//parm.Cluster_Two_Lines_X (0 . 25, 0. 750) ; 




031 




//parm.Cluster_Two_Lines_Y (0 . 25, 0. 150) ; 




032 




//Example 3 




033 




parm.Bound_Clust_X (0.5) ; 




034 




parm.Bound_Clust_Y (0.5) ; 




035 








036 




parm.Fill_del_xi_eta (del_xi, del_eta) ; 




037 








038 




Domain physical (xdim, ydim) ; 




039 




physical .Read_Bd ( ) ; 




040 




physical . Fill_del_xi_eta (del_xi, del_eta) ; 




041 




unsigned max_iter = 100; 




042 




double w = 1 . 90; 




043 




SORSOLVER {physical, parm, xdim, ydim) ; 




044 




//Reference or computational space 




045 




Domain ref (xdim, ydim) ; 




046 








047 




//Writing the mesh in the physical space ( GMV) 




048 




std: :ofstream outPhy ( "gmv_Physical . dat ",std: :ios: : out) 




04 9 




if ( ! (outPhy) ) std::cerr << "ERROR : UNABLE TO OPEN \" 


outPhy\' ' \n" ; 


050 




physical . GMV_Writer ( outPhy) ; 




051 




if (outPhy . is_open ( ) ) outPhy .close ( ) ; 




052 








053 




//writing mesh in the parameter space (GMV) 




054 




std : : of stream out Parm ( "gmv_Para . dat" , std : : ios : : out ) ; 




055 




if ( ! (outParm) ) std::cerr << "ERROR : UNABLE TO OPEN \' 


' outParm\' ' \n 


056 




parm . GMV_Writ er { out Parm) ; 




057 




if (out Parm . is_open ( ) ) out Parm .closed ; 




058 




parm.Matlab_Writer () ; 




059 








060 




return EXIT_SUCCESS; 




061 


} 






62 
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063 



2.2 domain.h 

001 #ifndef PARAMETER_SPACE 

002 #define P AR AME TE R_S P AC E 

004 #include<iost ream> 

005 #include<iomanip> 
00 6 #include<vector> 

007 #include<iterator> 

008 #include<f stream> 
00 9 #include<sstream> 
010 #include<map> 

012 # include "matrix . h" 

014 class Domain{ 

015 public: 



016 Domain ( ) ; 

017 Domain (unsigned int xdiml, unsigned ydiml) ; 

018 Domain (const Domain & org) ; 

019 unsigned int XDIM() const ; 
02 unsigned int YDIM ( ) const ; 

021 void Grid_Gen(); 

022 std: : vector<double> XCOORDS ( ) ; 

023 std: : vector<double> YCOORDSf); 

024 double Eriksson_l (double eta); 

025 

026 void Cluster_X_Near (double etaO); 

027 void Cluster_Y_Near (double etaO); 

028 

029 void Cluster_Two_Lines_X (double etal, double eta2); 

030 void Cluster_Two_Lines_Y (double etal, double eta2); 

031 

032 void Bound_Clust_X (double etal); 

033 void Bound_Clust_Y (double etal); 

034 

035 doubles XCOORD (unsigned int i , unsigned j); 

036 doubles YCOORD (unsigned int i , unsigned j); 

037 void Read_Bd ( ) ; 

038 void Matlab_Writer ( ) ; 

039 void GMV_Writer (std: : of stream & outFile) ; 

040 void Fill_del_xi_eta (double xi, double eta); 

041 //void Call_Grid_Adapter () ; 
042 

044 std: :vector<double> Pll (unsigned int i , unsigned int j); 

045 std: :vector<double> P22 (unsigned int i , unsigned int j); 

046 std: :vector<double> P12 (unsigned int i , unsigned int j); 

047 //++++++++++ 

048 Matrix MeshX () ; 

049 Matrix MeshY () ; 

051 double G22 (unsigned int i , unsigned int j); 

052 double Gil (unsigned int i , unsigned int j); 

053 double X_xi (unsigned int i , unsigned int j); 

054 double Y_xi (unsigned int i , unsigned int j); 

055 double X_eta (unsigned int i , unsigned int j); 

056 double Y_eta (unsigned int i , unsigned int j); 

057 double X_xieta (unsigned int i , unsigned int j) ; 

058 double Y_xieta (unsigned int i , unsigned int j); 
059 

60 private: 

061 double del_eta, del_xi; 

062 unsigned xdim, ydim; 

063 Matrix x, y; 

064 std: : vector<double> xcoords , ycoords ; 
65 }; 



066 #endif 
067 

068 



2.3 domain.cpp 

001 # include "domain . h" 
002 

003 #ifndef _FUNCTIONS_ 

004 # in elude " functions . h" 

005 #endif 
006 
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007 #include <cassert> 
008 

009 Domain :: Domain () { 

010 xdim = ; ydim = 0; 

011 } 

012 Domain :: Domain (unsigned int xdiml , unsigned int ydiml){ 

013 xdim = xdiml ; ydim = ydiml; 

014 } 

015 Domain :: Domain (const Domain s org) { 

016 xdim = org.XDIMO; 

017 ydim = org . YDIM ( ) ; 

018 Grid_Gen(); 

019 } 

020 unsigned int Domain :: XDIM ( ) const{ 

021 return xdim; 

022 } 

023 unsigned int Domain :: YDIM ( ) const{ 

024 return ydim; 

025 } 

026 void Domain: :Grid_Gen () { 

027 Matrix xt (xdim, ydim) , yt (xdim, ydim) ; 

028 assert(0 != xdim && != ydim); 

029 for (unsigned int j = ; j < ydim ; ++j){ 

030 for (unsigned int i = ; i < xdim ; ++i){ 

031 double t_x = double (i) /double (xdim-1 . 0) ; 

032 double t_y - double ( j ) /double (ydim-1 . ) ; 

033 xt(i,j) - t_x ; yt(i,j) - t_y; 

034 } 

035 } 

036 x - xt ; y = yt; 

037 } 

038 std: :vector<double> Domain :: XCOORDS () { 

039 xcoords . resize (xdim*ydim) ; ycoords . resize (xdim*ydim) ; 

040 for(int j = ; j < ydim ; ++j){ 

041 for (int i = ; i < xdim ; ++i){ 

042 int no = i+j*xdim; 

043 xcoords[no] = x(i,j); 

044 ycoords [no] = y(i,j); 

045 } 

046 } 

047 return xcoords; 

048 } 

049 std: :vector<double> Domain :: YCOORDS () { 

050 xcoords . resize (xdim*ydim) ; ycoords . resize (xdim*ydim) ; 

051 for(int j = ; j < ydim ; 

052 for (int i = ; i < xdim ; ++i){ 

053 int no = i+j*xdim; 

054 xcoords [no] = x(i,j); 

055 ycoords [no] = y(i,j); 

056 } 

057 } 

058 return ycoords; 

059 } 
060 

061 void Domain: :Bound_Clust_X (double etal){ 
062 

063 double alpha - 4.0; 

064 double h - 1.0; 

065 double h2 ■ 1.0; 

066 double hi - 0.0; 

067 

068 forfint j = ; j < ydim ; ++j){ 

069 forfint i = ; i < xdim ; ++i){ 

070 

071 if(x(i,j) <= etal && <= x(i,j)){ 

072 double eta = x(i,j); 

073 x(i,j) -(h2-hl)* et al* ( std :: exp (alpha*eta/etal ) -1 . )/( std :: exp ( alpha ) -1 . ) +hl ; 

074 } 

075 

076 if(x(i,j) >= etal ss x(i,j) <- 1.0){ 

077 double eta - x(i,j); 

078 x(i,j) - (h2-hl)*(1.0- (1.0-etal) * ( ( (std: : exp ( alpha* ( 1 . 0-eta) / (1.0-etal) ) ) -1 . 0) / (std: : exp (alpha) -1 .0) ) ) ; 

079 } 

080 

081 } 

082 } 
083 

084 } 

085 

086 void Domain: :Bound_Clust_Y (double etal){ 



087 

088 double alpha - 4.0; 

8 9 double h - 1.0; 

090 double h2 - 1.0; 

091 double hi - 0.0; 

092 

093 for(int j = ; j < ydim ; ++j){ 



094 forfint i = ; i < xdim ; ++i){ 
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095 

096 if(y(i,j) <- etal ss <- y(i,j)){ 

097 double eta - y(i,j); 

098 y(i,j) -(h2-hl)* et al* ( std :: exp (alpha*eta/etal ) -1 . )/( std :: exp ( alpha ) -1 . ) +hl ; 

099 } 

100 

101 if(y(i,j) >= etal ss y(i,j) <= 1.0){ 

102 double eta - y(i,j); 

103 y(i,j) - (h2-hl) * (1.0- (1.0-etal) * ( ( (std: : exp ( alpha* ( 1 . 0-eta) / (1.0-etal) ) ) -1 . 0) / (std: : exp (alpha) -1 .0) ) ) ; 

104 } 

105 

106 } 

107 } 

108 



109 } 

110 

111 //———-- 

112 double Domain :: Eriksson_l (double eta) { 

113 double h - 1.0; 

114 double alpha - 3.0; 

115 

116 return h* (( std :: exp ( alpha*eta) -1 . )/( std :: exp (alpha) -1 . )) ; 

117 
118 
119 } 
120 

121 void Domain: : Cluster_Two_Lines_X (double etal, double eta2){ 



122 

123 double alpha = 5.0; 

124 double h - 1.0; 

125 double etaO = ( et al+et a2 ) *0 . 5 ; 
126 

127 for(int j = ; j < ydim ; ++j){ 

128 for(int i = ; i < xdim ; ++i){ 
129 

130 if(x(i,j) <- etal ss <- x(i,j)){ 

131 double eta - x(i,j); 

132 x(i,j) - etal* (h-Eriksson_l (1-eta/etal) ) ; 

133 } 

134 

135 if(x(i,j) >= etal ss x(i,j) <= etaO ){ 

136 double eta - x(i,j); 

137 x(i,j) - h*etal+ ( et aO-et al ) *Eriksson_l ( (eta-etal ) / (etaO-etal ) ) ; 

138 

139 } 
140 

141 if(x(i,j) >= etaO ss x(i,j) <= eta2){ 

142 double eta - x(i,j); 

143 x(i,j) - h*etaO + (eta2-eta0) * (h-Eriksson_l ( (eta2-eta) / (eta2-eta0) ) ) ; 

144 } 

145 

146 if(x(i,j) >= eta2 ss x(i,j) <- 1.0){ 

147 double eta - x(i,j); 

148 x(i,j) - h*eta2 + ( 1 . 0-et a2 ) *Eriksson_l ( (eta-eta2 ) / ( 1 . 0-eta2 ) ) ; 

149 } 

150 

151 } 

152 } 

153 



154 } 

155 

156 void Domain :: Cluster_Two_Lines_Y (double etal, double eta2){ 



157 

158 double alpha = 5.0; 

15 9 double h - 1.0; 

160 double etaO - (etal+eta2 ) *0 . 5; 

161 

162 for(int j = ; j < ydim ; ++j){ 

163 for(int i = ; i < xdim ; ++i){ 
164 

165 if(y(i,j) <= etal ss <- y(i,j)){ 

166 double eta = y(i,j); 

167 y(i,j) - etal* (h-Eriksson_l (1-eta/etal) ) ; 

168 } 
169 

170 if(y(i,j) >- etal ss y(i,j) <- etaO ){ 

171 double eta - y(i,j); 

172 y(i,j) - h*etal+ (etaO-etal) *Eriksson_l ( (eta-etal) / (etaO-etal) ) ; 

173 

174 } 

175 

176 if(y(i,j) >= etaO ss y(i,j) <= eta2){ 

177 double eta = y(i,j); 

178 y(i,j) - h*etaO + (et a2-et aO ) * (h-Eriksson_l ( (eta2-eta) / (eta2-eta0 ) ) ) ; 

179 } 

180 

181 if(y(i,j) >= eta2 ss y(i,j) <= 1.0){ 

182 double eta - y(i,j); 
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183 y(i,j) - h*eta2 + ( 1 . O-et a2 ) *Eriksson_l ( (eta-eta2 ) / ( 1 . 0-eta2 ) ) ; 

184 } 

185 

186 } 

187 } 
188 

189 } 

190 //============= 

191 void Domain: :Cluster_X_Near (double eta0){ 



192 double alpha - 3.0; 

193 for(int j = ; j < ydim ; ++j){ 

194 for (int i = ; i < xdim ; ++i){ 

195 if (x (i, j) < etaO) { 

196 double eta = x(i,j); 

197 x(i,j) - (double) etaO * (exp(alpha) - exp(alpha * (double) (1 - eta / etaO))) / (exp(alpha) - O.lel); 

198 } 

199 if (x (i, j) > etaO) { 

200 double eta - x(i,j); 

201 x(i,j) - (double) etaO + (double) (1 - etaO) * (exp ( (double) (alpha * (eta - etaO) / (1 - etaO))) - O.lel) 

202 / (exp ( (double) alpha) - O.lel); 

203 } 

204 } 

205 } 
206 



207 } 
208 

209 void Domain: :Cluster_Y_Near (double eta0){ 



210 double alpha - 3.0; 

211 forfint j = ; j < ydim ; ++j){ 

212 forfint i = ; i < xdim ; ++i){ 

213 if (y (i, j) < etaO) { 

214 double eta = y(i,j); 

215 y(i,j) - (double) etaO * (exp(alpha) - exp(alpha * (double) (1 - eta / etaO))) / (exp(alpha) - O.lel); 

216 } 

217 if (y (i, j) > etaO) { 

218 double eta - y(i,j); 

219 y(i,j) - (double) etaO + (double) (1 - etaO) * (exp ( (double) (alpha * (eta - etaO) / (1 - etaO))) - O.lel) 

220 / (exp ( (double) alpha) - O.lel); 

221 } 

222 } 

223 } 

224 } 

225 doubles Domain :: XCOORD (unsigned int i , unsigned int j) { 

226 if (i >- xdim I I j >- ydim || i<0 || j < 0) { 

227 std: icerr << "In XCOORDS (..,.. ) dim mismatch\n"; 

228 } 

229 return x (i, j) ; 

230 } 

231 doubles Domain :: YCOORD (unsigned int i , unsigned int j){ 

232 if (i >- xdim I I j >- ydim || i<0 || j < 0) { 

233 std: icerr << "In YCOORDS (..,.. ) dim mismatch\n"; 

234 } 

235 return y (i, j) ; 

236 } 

237 void Domain: :Read_Bd() { 

238 //Read the boundary of the physical domain 

239 Matrix xt (xdim, ydim) , yt (xdim, ydim) ; 

240 forfint j - ; j < ydim ; ++j){ 

241 for (int i =0 ; i < xdim ; ++i){ 

242 if (0 — ill xdim-1 ~ i I I — j I I ydim-1 == j ) { 

243 double zetal = double ( i ) /double (xdim-1 . ) ; 

244 double etal = double ( j ) /double (ydim-1 . ) ; 

245 xt(i,j) = zetal; yt(i,j) = etal; 

246 

247 xt(i,j) - XYoirole (zetal, etal, 1) ; 

248 yt(i,j) - XYoirole (zetal, etal, 2) ; 
249 

250 } 

251 } 

252 } 

253 //Create grid by the TFI 

254 for(int j = 1 ; j < ydim-1 ; ++j){ 

255 forfint i = 1 ; i < xdim-1 ; ++i){ 

256 double zetal = double ( i ) /double (xdim-1 . ) ; 

257 double etal - double ( j ) /double (ydim-1 . ) ; 

258 xt(i,j) - (1 . 0-zetal) *xt (0, j) +zetal*xt (xdim-1 , j ) + ( 1 . 0-etal ) *xt ( i, ) +etal*xt ( i , ydim-1 ) - 

259 ( (1. 0-zetal) * (1 ■ 0-etal) *xt (0, 0) + ( zetal )*( 1 . 0-etal ) *xt (xdim-1 , ) 

260 + (zetal) * (etal) *xt (xdim-1, ydim-1) + (1 . 0-zetal ) *etal*xt (0, ydim-1 )) ; 

261 yt(i,j) - (1-zetal) *yt (0, j) +zetal*yt (xdim-1, j) + ( 1 . 0-etal ) *yt ( i , ) +etal*yt ( i, ydim-1 ) - 

262 ( (1-zetal) * (1-etal) *yt (0, 0) + ( zetal )*( 1-etal ) *yt (xdim-1 , ) 

263 + (zetal) * (etal) *yt (xdim-1, ydim-1) + ( 1-zetal ) *etal*yt ( , ydim-1 )) ; 

264 } 

265 } 

2 66 x - xt ; y = yt ; 



267 } 
2 68 //============= 

269 void Domain : : Mat lab_Writ er () { 

270 
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271 std: : vector<double> xl = XCOORDSf); 

272 std: : vector<double> yl = YCOORDSf); 

273 std: : vector<double> : : const_iterator viter; 

27 4 std: : of stream out file ( "matlab_out .in", std: :ios: : out ) ; 

275 if(loutfile) std::cerr << "Unable to open the matlab outfile\n"; 

276 outfile << "clear;\n"; 

277 outfile << "holdon=ishold; \n" ; 
278 

279 for(int j = ; j < ydim ; 

280 for(int i = ; i < xdim ; ++i){ 
2 81 int no = i + j *xdim; 

282 outfile << "xl(" << i+1 << "," << j+1 << ")=" << xl [no] << "; 

283 << "yl(" << i+1 << "," << j+1 << ")=" << yl[no] << ";" << std::endl; 

284 } 

285 } 
286 

287 outfile << "m = " << xdim << std: :endl; 

288 outfile << "n = " << ydim << std: :endl; 
289 

290 outfile << "plot (xl ( 1 , : ) , y 1 ( 1 , : ) , ' r ' ) ; hold on" << std::endl; 

2 91 outfile << "plot (xl (m, : ) , yl (m, : ) , ' r' ) ; " << std:: endl; 

292 outfile << "plot (xl ( : , 1) , yl ( : , 1) , ' r' ) ; " << std::endl; 

293 outfile << "plot (xl ( : , n) , yl ( : , n) , ' r' ) ; " << std::endl; 

294 

295 outfile << "% Plot internal grid lines\n" ; 

296 outfile << "for i=2:m-l, plot (xl ( i , : ) , yl ( i , : ) , ' b' ) ; end\n"; 

297 outfile << "for j=2:n-l, plot (xl ( : , j ) , yl ( : , j ) , ' b' ) ; end\n"; 

298 

299 outfile << "if (~holdon), hold off, end" << std: : endl; 
300 

301 outfile << "axis off;\n"; 

302 

3 03 outfile .closed ; 
304 } 

305 

306 void Domain :: GMV_Writer (std: : of stream & outFile) { 
307 

308 std: :vector<double> xcoords = XCOORDS () ; 

309 std: :vector<double> ycoords = YCOORDSf); 
310 

311 outFile << "gmvinput ascii\n " ; 

312 outFile << "nodes " << xdim* ydim << std : : endl ; 
313 

314 for (int j = ; j < ydim ; ++j){ 

315 for (int i = ; i < xdim ; + + i){ 

316 int no = i + j*xdim; 

317 outFile << xcoords [no] << " "; 
318 

319 } 

320 } 

321 outFile << std::endl << std::endl; 

322 //writing y coord 

323 for (int j = ; j < ydim ; ++j){ 
32 4 for (int i = ; i < xdim ; + + i) { 
32 5 int no = i + j*xdim; 

326 outFile << ycoords [no] << " " ; 

327 

328 } 

329 } 

330 outFile << std::endl << std::endl; 

331 //forming cells 

332 outFile << "cells " << (xdim-1 ) * ( ydirn-1 ) << std::endl; 

333 for (int j = ; j < (ydim-1) ; ++j){ 

334 for (int i = ; i < (xdim-1) ; + + i) { 

335 int no = ( i + (xdirn) ) +1; 

336 int nol = ( i + ( j+1 )* (xdim) +1 ) ; 

337 outFile << "quad 4 " << std::endl; 

338 outFile << no << " " << no+1 << " " 

339 << nol+1 << " " << nol << std: : endl; 

340 } 

341 } 

342 outFile << std: : endl; 
343 

344 outFile << std::endl << "endgmv\n"; 

345 outFile. close () ; 

346 } 

348 Matrix Domain :: MeshX () { 

34 9 return x; 

350 } 

351 Matrix Domain :: MeshY () { 
352 

3 53 return y; 

354 } 

356 double Domain :: G22 (unsigned int i , unsigned int j){ 

3 57 double xl ; double x2 ; double y 1 ; double y2 ; 

358 xl = x(i,j-l); x2 = x(i,j + l) ; 
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359 yl - y(i,j-l); y2 - y(i,j+l); 

360 double g22 - st d : : pow ( (x2-xl ) / ( 2 . 0*del_et a ) , 2 ) + 

361 std: :pow( (y2-yl) / ( 2 . tdel_et a ) , 2) ; 
3 62 return g22; 

363 } 

364 double Domain :: Gil (unsigned int i , unsigned int j){ 

365 double xl = x(i-l,j); double x2 = x(i+l,j); 

366 double yl - y(i-l,j); double y2 - y(i+l,j); 

367 double gll - std : : pow ( (x2-xl ) / (2 . 0*del_xi) , 2 ) + 
3 68 std: :pow( (y2-yl) / (2.0*del_xi) , 2) ; 

369 return gll; 

370 } 

371 double Domain :: X_xi (unsigned int i , unsigned int j){ 
3 72 double x_xi; 

373 x_xi - (x (i + 1, j) -x (i-1, j) ) / (2 . 0*del_xi) ; 

37 4 return x_xi; 

375 } 

376 double Domain :: X_et a (unsigned int i , unsigned int j){ 

377 double x_eta; 

378 x_eta - (x ( i, j+1 ) -x ( i, j-1 ) ) / ( 2 . 0*del_et a ) ; 
37 9 return x_eta; 

380 } 

381 double Domain :: Y_xi (unsigned int i , unsigned int j ) { 
3 82 double y_xi; 

383 y_xi - (y (i+1, j) -y (i-1, j) ) / (2 .0*del_xi) ; 

3 8 4 return y_xi ; 

385 } 

386 double Domain :: Y_eta (unsigned int i , unsigned int j ) { 

387 double y_eta; 

388 y_eta - ( y ( i, j + 1 ) -y ( i, j-1 ) ) / ( 2 . 0*del_et a ) ; 

389 return y_eta; 

390 } 

391 double Domain :: X_xiet a (unsigned int i , unsigned int j){ 

3 92 return (x (i + 1 , j + 1 ) +x (i-1 , j-1 ) -x (i-1 , j + 1 ) -x (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_et a ) ; 

393 } 

394 double Domain :: Y_xieta (unsigned int i , unsigned int j){ 

3 95 return (y (i + 1 , j + 1 ) +y (i-1 , j-1 ) -y (i-1 , j + 1 ) -y (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_et a ) ; 

396 } 

397 void Domain :: Fill_del_xi_eta (double xi, double eta) { 

398 del_xi = xi; del_eta = eta; 

399 } 

400 std: :vector<double> Domain :: P 1 1 (unsigned int i , unsigned int j){ 

401 //x-t coordinate 

402 //compute the jacobian at the point 

403 double s_xi - (x ( i + 1 , j ) -x ( i-1 , j ) ) / ( 2 . 0*del_xi ) ; 

404 double s_eta - (x (i, j + 1 ) -x (i, j-1 ) ) / (2 . 0*del_eta) ; 

405 double t_xi - (y (i+1 , j ) -y (i-1 , j ) ) / (2 . 0*del_xi) ; 

406 double t_eta - (y (i, j+1) -y (i, j-1) ) / (2 . 0*del_eta) ; 

407 double det = s_xi*t_eta-t_xi*s_eta; 

408 double TI_11 = t_eta/det; double TI_12 = -t_xi/det; 

409 double TI_21 = -s_eta/det; double TI_22 = s_xi/det; 

410 double s_xixi = (x ( i+ 1 , j ) -2 . 0*x ( i , j ) +x ( i-1 , j ) ) / (del_xi*del_xi ) ; 

411 double s_etaeta - (x ( i , j +1 ) -2 . 0*x ( i, j ) +x ( i, j-1 ) ) / ( del_et a*del_et a ) ; 

412 double s_xieta = (x (i+1 , j+1 ) +x (i-1 , j-1 ) -x (i-1 , j+1 ) -x (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

413 double t_xieta = (y (i+1 , j + 1 ) +y (i-1 , j-1 ) -y (i-1 , j+1 ) -y (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

414 double t_xixi = (y (i+1 , j ) -2 . 0*y (i, j ) +y (i-1 , j ) ) / (del_xi*del_xi ) ; 

415 double t_etaeta - (y ( i , j + 1 ) -2 . 0*y ( i, j ) +y ( i, j-1 ) ) / (del_et a*del_et a ) ; 

416 std: :vector<double> Pll(2); 
417 

418 P11[0] = - (s_xixi*TI_ll+t_xixi*TI_12) ; 

419 Pll[l] = - (s_xixi*TI_21+t_xixi*TI_22) ; 

420 

421 return Pll; 

422 } 

423 std: :vector<double> Domain :: P22 (unsigned int i , unsigned int j){ 

424 //x-t coordinate 

425 //compute the jacobian at the point 

426 double s_xi - (x (i+1 , j ) -x (i-1 , j ) ) / (2 . 0*del_xi) ; 

427 double s_eta - (x (i, j+1 ) -x (i, j-1 ) ) / (2 . 0*del_eta) ; 

428 double t_xi - (y (i+1 , j ) -y (i-1 , j ) ) / (2 . 0*del_xi) ; 

429 double t_eta - (y (i, j + 1 ) -y (i, j-1 ) ) / (2 . 0*del_eta) ; 

430 double det = s_xi*t_eta-t_xi*s_eta; 

431 double TI_11 = t_eta/det; double TI_12 = -t_xi/det; 

432 double TI_21 = -s_eta/det; double TI_22 = s_xi/det; 

433 double s_xixi = (x ( i+ 1 , j ) -2 . 0*x ( i , j ) +x ( i-1 , j ) ) / (del_xi*del_xi ) ; 

434 double s_etaeta - (x ( i , j +1 ) -2 . 0*x ( i, j ) +x ( i, j-1 ) ) / ( del_et a*del_et a ) ; 

435 double s_xieta - (x (i+1 , j+1 ) +x (i-1 , j-1 ) -x (i-1 , j+1 ) -x (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

436 double t_xieta - (y (i+1 , j + 1 ) +y (i-1 , j-1 ) -y (i-1 , j+1 ) -y (i+1 , j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

437 double t_xixi - (y ( i + 1 , j ) -2 . 0*y ( i, j ) +y ( i-1 , j ) ) / ( del_xi*del_xi ) ; 

438 double t_etaeta - (y ( i , j +1 ) -2 . 0*y ( i, j ) +y ( i, j-1 ) ) / (del_et a*del_et a ) ; 

439 std: :vector<double> P22(2); 

440 P22[0] - - (s_etaeta*TI_ll+t_etaeta*TI_12) ; 

441 P22[l] - - (s_etaeta*TI_21+t_etaeta*TI_22) ; 

442 return P22; 

443 } 

444 std: :vector<double> Domain :: P 12 (unsigned int i , unsigned int j){ 

445 //x-t coordinate 

446 //compute the jacobian at the point 
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447 double s_xi = (x (i+1 , j ) -x (i-1 , j ) ) / (2 . 0*del_xi ) ; 

448 double s_eta - (x (i, j+1 ) -x (i, j-1 ) ) / (2 . 0*del_eta) ; 

449 double t_xi - (y (i+1 , j ) -y (i-1 , j ) ) / (2 . 0*del_xi ) ; 

450 double t_eta - (y (i, j+1) -y (i, j-1) ) / (2 . 0*del_eta) ; 

451 double det = s_xi*t_eta-t_xi*s_eta; 

452 double TI_11 = t_eta/det; double TI_12 = -t_xi/det; 

453 double TI_21 = -s_eta/det; double TI_22 = s_xi/det; 

454 double s_xixi = (x ( i + 1 , j ) -2 . 0*x ( i , j ) +x ( i-1 , j ) ) / (del_xi*del_xi ) ; 

455 double s_etaeta - (x (i , j+1 ) -2 . 0*x (i, j ) +x (i, j-1 ) ) / (del_eta*del_eta) ; 

456 double s_xieta - (x (i+1 , j + 1 ) +x (i-1 , j-1 ) -x (i-1 , j+1 ) -x (i+1, j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

457 double t_xieta = (y (i + 1 , j + 1 ) +y (i-1 , j-1 ) -y (i-1, j+1 ) -y (i + 1 , j-1 ) ) / ( 4 . 0*del_xi*del_eta) ; 

458 double t_xixi = (y (i+1 , j ) -2 . 0*y (i, j ) +y (i-1 , j ) ) / (del_xi*del_xi) ; 

459 double t_etaeta - (y (i, j+1 ) -2 . 0*y (i, j ) +y (i, j-1 ) ) / (del_eta*del_eta) ; 

460 std: :vector<double> P12(2); 

461 P12[0] = - (s_xieta*TI_ll+t_xieta*TI_12) ; 

462 P12[l] = - (s_xieta*TI_21+t_xieta*TI_22) ; 

463 return P12; 



464 } 
465 

466 



2.4 matrix.cpp 

001 #ifndef MATRIX_H 

002 #define MATRIX_H 

003 class Matrix{ 

004 public: 



005 unsigned int nx,ny; 

006 std: :vector<std: : vector<double> > Elements; 

007 Matrix (){ 

008 } 
009 

010 Matrix (unsigned int nxl, unsigned int nyl){ 

011 nx = nxl ; ny = nyl; 

012 //number of rows 

013 Elements . resize (nx) ; 

014 //fill each rows 

015 for (unsigned int i = ; i < nx ; ++i) 

016 Elements [i] . resize (ny, 99 . 0) ; 

017 } 

018 

019 void Clear () { 

020 for (unsigned int i = ; i < nx ; ++i) 

021 Elements [i] . clear () ; 

022 } 

023 

024 doubles operator ( ) (unsigned int ix, unsigned int iy){ 

025 return Element s [ ix] [ iy] ; 

026 } 



027 }; 

028 tendif 
029 

030 



2.5 sor_solver.cpp 

001 #ifndef SOR_SOLVER 

002 #define SOR_SOLVER 
003 

004 #include <vector> 
005 

00 6 fin elude "domain . h" 
007 # include "matrix . h" 
008 

009 double Mesh_Residual (Matrix x, Matrix y , 



010 Matrix x_old, Matrix y_old, 

011 unsigned int xdim, unsigned int ydim) { 

012 double resid = ; 

013 for (int j = ; j < ydim ; ++j){ 

014 for (int i = ; i < xdim ; ++i) { 

015 double x_resd = (x ( i , j ) -x_old ( i, j ) ) ; 

016 double y_resd = (y ( i , j ) -y_old ( i, j ) ) ; 

017 resid += (x_resd*x_resd+y_resd*y_resd) ; 

018 } 

019 } 

020 return std: : sqrt (resid) ; 



021 } 
022 

023 bool S0RS0LVER( Domains physical, Domain & parm , unsigned int xdiml , unsigned int ydiml) { 
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024 bool grid_dist = true; 

025 bool run_ellip = true; 

026 

027 double tolerance = 1.0e-4; 

028 unsigned max_iter = 100; 

029 double w = 1.90; 

030 double residual = 10.0; 

031 unsigned int iter = 0; 
032 

033 unsigned int xdim =xdiml ; unsigned int ydim = ydiml ; 

03 4 Matrix x_old (xdim, ydim) , y_old (xdim, ydim) ; 

035 

036 double del_xi = 1 . 0/double (xdim-1 . 0) ; 

037 double del_eta = 1 . 0/double (ydim-1 . 0) ; 

03 8 

03 9 std: :ofst ream f out ( "gauss . dat " , std: : ios : : out) ; 

040 if ( ! (fout .is_open () ) ) std::cerr << "ERROR : UNABLE TO OPEN THE FILE \" gauss . dat\" \n" ; 
041 

042 if (run_ellip) { 

043 while (iter < max_iter & residual > tolerance) { 

044 iter++; 

4 5 x_old = physical .MeshX ( ) ; y_old = physical . MeshY ( ) ; 

046 for (unsigned int j = 1 ; j < ydim-1 ; ++j){ 

047 for (unsigned int i = 1 ; i < xdim-1 ; ++i) { 

048 

049 double g22 = physical . G22 ( i , j ) ; 

050 double gll = physical . Gil ( i , j ) ; 

051 double x_xi = physical . X_xi ( i, j ) ; 

52 double x_et a = physical . X_et a ( i, j ) ; 

53 double y_xi = physical . Y_xi ( i, j ) ; 

054 double y_eta = physical . Y_et a ( i, j ) ; 

055 double x_xieta = physical . X_xiet a ( i, j ) ; 

056 double y_xieta = physical . Y_xiet a ( i, j ) ; 

057 double gl2 = x_xi*x_eta+y_xi*y_eta; 

058 double g = std : : pow (x_xi*y_eta-y_xi*x_eta, 2 ) ; 
059 

060 std: :vector<double> Pll , P22,P12; 
061 

062 if (grid_dist) { 

063 Pll = parm.Pll (i, j) ; 

064 P22 = parm.P22 (i, j) ; 

065 P12 = parm.P12 (i, j) ; 

066 }else{ 

067 Pll .push_back (0) ; P 1 1 . push_back ( ) ; 
68 P22 .push_back (0) ; P22 . push_back ( ) ; 
69 P12 .push_back (0) ; P 12 . push_back ( ) ; 
070 } 

071 

072 double tmpx, tmpy ; 
073 

074 tmpx = (g22*Pll [0] -2 . 0*gl2*P12 [0] +gll*P22 [0] ) *x_xi+ 

07 5 (g22*Pll [1] -2 . 0*gl2*P12 [ 1 ] +gl 1*P22 [ 1 ] ) *x_eta; 

076 

077 tmpy = (g22*Pll [ ] -2 . 0*gl2*P12 [ ] +gll*P22 [ ] ) *y_xi+ 

078 (g22*Pll [1] -2 . 0*gl2*P12 [ 1 ] +gl 1*P22 [ 1 ] ) *y_eta; 

079 

080 double lhsx = 2 . 0* (g22/ ( del_xi*del_xi ) +gll/ (del_eta*del_eta) ) ; 

081 double rhsx = g22* (physical . XCOORD ( i + 1 , j ) +physical . XCOORD ( i-1 , j ))/( del_xi*del_xi ) 
82 +gll* (physical .XCOORD (i, j + 1 ) +physical . XCOORD ( i, j-1) ) / ( del_et a*del_et a ) 

083 -2 . 0*gl2*x_xieta + tmpx; 
084 

085 physical .XCOORD (i, j) = physical . XCOORD ( i, j ) + w* (rhsx/ lhsx-physical . XCOORD ( i , j ) ) ; 
086 

087 double lhsy = lhsx; 

088 double rhsy = g22* (physical . YCOORD ( i + 1 , j ) +physical . YCOORD ( i-1 , j ))/( del_xi*del_xi ) 

089 +gll* (physical .YCOORD (i, j + 1 ) +physical . YCOORD ( i, j-1) ) / ( del_et a*del_et a ) 

090 -2 . 0*gl2*y_xieta+ tmpy; 
091 

092 physical .YCOORD (i, j) = physical . YCOORD ( i, j ) + w* (rhsy/ lhsy-physical . YCOORD ( i , j )) ; 

093 } 

094 } 
095 

96 double mesh_resid = Mesh_Residual (physical .MeshX ( ) , physical . MeshY ( ) , x_old, y_old, xdim, ydim) 

097 / ( (xdim-2) * (ydim-2) ) ; 

098 std::cout << "Iteration = " << iter << ", Residual = " 

099 << mesh_resid << std::endl; 

100 

101 fout << iter << " " << mesh_resid << std::endl; 

102 residual = mesh_resid; 

103 } 

104 } 
105 

10 6 return true; 



107 } 

108 #endif 

109 

110 
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2.6 functions.h 



001 #ifndef _FUNCTIONS_ 

002 #define _FUNCTIONS_ 
003 

004 #include <vector> 

005 #include <cmath> 
00 6 #include<iomanip> 

07 # include <iost ream> 
008 

9 double XYcircle (double x, double y , unsigned int x_or_y ) { 



010 double r = 1.0; 

011 double theta = 0.0; 

012 double Pi = 4 . 0*atan (1 . 0) ; 

013 if (0==y) { 

014 theta = Pi/2.0*x; 

015 if(x_or_y == 1) 

016 return r*cos (theta) ; 

017 else 

017 return r*sin (theta) ; 

018 } 

019 if (l==x) { 

020 theta = Pi/2+Pi/2*y; 

021 if(x_or_y == 1) 

022 return r*cos (theta) ; 

023 else 

023 return r*sin (theta) ; 

024 } 

025 if(l==y){ 

02 6 theta=Pi+Pi/2* (1 . 0-x) ; 

027 if(x_or_y == 1) 

028 return r*cos (theta) ; 

029 else 

029 return r*sin (theta) ; 

030 } 

031 if (0==x) { 

032 theta = 3 . 0*Pi/2 . 0+Pi/2 . 0* (1 . 0-y) ; 

033 if(x_or_y == 1) 

034 return r*cos (theta) ; 

035 else 

035 return r*sin (theta) ; 

036 } 



037 } 

038 #endif 
039 

040 



2.7 makefile 

001 CXX - g++ 

002 SRC = main.cpp\ 

003 domain . cpp 

004 OBJ = main.o\ 

005 domain . o 

006 ell_mesh: $ (OBJ) 

007 $ (CXX) -Wall -03 -o ellmesh $ (OBJ) $ (LIB) 

008 clean: 

009 rm -f ellmesh 

010 rm -f *.o 

on 

012 



3 Numerical Examples 



3.1 Example 1 

In this example, we cluster the grids along the centre lines of a circular physical domain. Figure 
|4]shows the grid in the parameter space for concentrating grids at the centre lines of the physical 
space. Grid density in the physical space is determined by the grid density in the parameter 
space. Figure |5] shows the converged grid in the physical space. For generating this grid the 
lines 027 and 028 in the subsection l2.1l are used. 
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Figure 4: Grid in the parameter space. 



Figure 5: Adapted grid by elliptic system. 




Figure 6: Grid in the parameter space. Figure 7: Adapted grid by elliptic system. 

3.2 Example 2 

See the Figure|6]for grids in the parameter space and the Figure|7]for the converged grids in the 
physical space. For generating the grids the lines 030 and 031 of the subsection l2.1l are used. 



3.3 Example 3 

In this example, we are interested in concentrating grids at the boundary of the physical space. 
Figure [U shows the grid in the parameter space for concentrating grids at the boundary of the 
physical domain. The converged grids in the physical space is shown in the Figure For 
generating the grids the lines 033 and 034 of the subsection l2.1l are used. 
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Figure 8: Grid in the parameter space. Figure 9: Adapted grid by elliptic system. 

4 Conclusions 

An elliptic system for generating adaptive quadrilateral meshes in curved domains has been 
presented. A C ++ implementation of the presented technique is also given. Three examples are 
reported for demonstrating the effectiveness of the technique and the implementation. Since, 
the quadrilateral meshes are very extensively used for numerical simulations and grid adaptation 
is required for capturing many important phenomenon such as the boundary layers. Thus, this 
software package can be a very useful tool. 

Our package is freely available at www.mi.uib.no/~sanjay. The authors want to 
mention that they donot know of any freely available software package that can generate adap- 
tive quadrilateral meshes. 
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